Quasiparticle dynamics by effective π\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi $$\end{document}-field distortion

Modeling dynamical processes of quasiparticles in low dimensional π\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi $$\end{document}-conjugated systems is challenging due to electron-phonon coupling. We show that this interaction leads to linear potential energy terms in the lattice Lagrangian similar to a local “gravitational” field. The presence of quasiparticles deforms this field in a way analogous to a low-dimension solution of general relativity. Our calculations with analytical expressions for effective π\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi $$\end{document}-fields yield the correct band structure and deliver proper time evolution of the quasiparticle’s properties. Furthermore, we report a sharp reduction in the dynamics computational time up to two orders of magnitude, a result that has major simulation implications.

www.nature.com/scientificreports/ other phenomena reveals the equivalence between the methodologies. Moreover, the method is also notably cost-efficient, producing time reductions of two orders of magnitude, indicating its use for larger systems in which the conventional procedure becomes impracticable.

Methods
We present the SSH model in detail to point out precisely what can be done in a different way. Cis-polyacetylene is considered to keep the notation as concise as possible. Besides, its structure is a common backbone of conducting polymers, and the results are extended easily to 2-D π-systems as graphene nanoribbons. In this formalism, the model Hamiltonian is the sum of electronic ( H tb ) and lattice ( H latt ) terms. The first contribution is a tightbinding-like operator given by 30 C † n+1,s is the creation operator of a π-electron at site n + 1 with spin s. t n+1,n is the hopping integral associated with sites n + 1 and n. The SSH model restricts the physical description to only low-energy excitations 14 . As a result, the σ-bond fluctuates close to its undisturbed length. The hopping integral is assumed to be dependent on the inter-site distance. Expanding t n+1,n over this length gives where t n ≡ (1 + (−1) n δ 0 )t 0 , with t 0 the hopping integral in the absence of distortions, δ 0 a Brazovskii-Kirova symmetry breaking term, and u n the displacement of site n with respect to its equilibrium position. The dependence of the hopping integral with the position is mediated by α , the electron-phonon coupling constant.
In addition, the small fluctuation of the σ-bonds allows to adopt the harmonic approximation and treat the lattice as a sum of coupled harmonic oscillators, Here, K stands for the elastic constant, M is the site's mass, and P i is the i-th conjugated momentum.
It is not possible to diagonalize the electronic Hamiltonian directly since it depends on the lattice configuration. Instead, one must simultaneously solve for both the electronic and lattice in a self-consistent manner 10 .
One begins by finding stationary solutions for the lattice and π-electrons in an iterative self-consistent way until the convergence criterion is matched. Once the solution converges, it is possible to evolve it in time by including an external electric field. It can be done through a Peierls substitution on the Hamiltonian 25 . Ultimately, this operation only incorporates a time-dependent phase to the hopping integral in Eq. 2, which now reads γ ≡ ea/( c) , where e stands for the elementary charge, a the lattice parameter, c the speed of light, and A(t) is the vector potential. The electric field depends only on time, E(t) = −(1/c)Ȧ(t).
Let |ψ k,s (t)〉 be the electronic solution at a given time t, and U(dt) ≡ e −iH(t)dt/ the time-evolution operator. Then, the state, after a time interval dt, becomes Expanding it over the basis of eigenstates, {|φ ls (t)�} , of the electronic Hamiltonian at a given time t leads to, or, As for the lattice, the expected value of the Lagrangian, L , obeys the Euler-Lagrange equations, L is calculated with a Slater-determinant state composed with all occupied electronic states, ψ k,s (n, t) . This equation leads to with, (2) t n+1,n = t n + α(u n+1 − u n ), (4) t n+1,n = e −iγ A(t) (t n + α(u n+1 − u n )).
(9) Mÿ n = K(y n+1 − 2y n + y n−1 ) +α(B n+1 − 2B n + B n−1 + c.c.), www.nature.com/scientificreports/ and, y n = u n+1 − u n . The prime sign means a summation over only occupied states. The time evolution for the model is determined by equations 7 and 9. When time advances one step, the electronic and lattice configurations change. Because of that, the set {φ ls (t ′ )} is not the same as before. After each time step, the modeling requires the re-diagonalization of the electronic Hamiltonian. Alternatively, one could use high-order explicit time integration methods to evolve the Schrödinger equation. However, these approaches have a significant impact on the algorithm's computational cost, meaning limitations on system size and simulation time. Equation 9 is the crucible of π-systems description. The B n terms depend on one-electron wavefunctions that are time demanding to find. But, B n can be calculated analytically for the stationary dimerized configuration 31 . In this case, the results render Eq. 9 as equivalent to a series of connected spring-masses subjected to an alternating gravitational field, since B n+1 = −B n . Which can also be deduced by symmetry arguments. Obviously, this picture changes as time and quasiparticles get into consideration. Nevertheless, the change is well-behaved and predictable. It might be more convenient to think in terms of a π-electrons force field defined by, Solving numerically in tandem Eqs. 7 and 9 for several systems 5,6,10 , it can be verified that far from quasiparticles and defects, n still behaves as an alternating effective "gravitational" field. It is also found that this effective field distorts in association with polarons according to the pattern sech 2 (x − vt) . It should be pointed out that the square hyperbolic secant dependence could be foreseen from the analytical solutions of the TLM continuum model 32 .
In what follows, we present the simulations of polaron dynamics with the usual approach and the proposed effective field to highlight the equivalence and difference between the two methods.
We set the parameters according to the usual values assigned to conjugated polymers. α = 4.1 eV/Å, t 0 = 2.5 eV, δ 0 = 0.05, K = 21 eV/Å 2 and a = 1.41 Å 19, 20 . We simulate a polaron state in the presence of an electric field with a strength of 2.4×10 5 V/cm. Here, E(t) will continuously grow from zero until t = 50 fs, when it reaches its maximum. Right after, the field drops to 0, remaining off until the end of the simulation. Throughout this work, we simulate the polymer with a system of 120 sites with periodic boundary conditions. The sites' masses equals 13 a.u.. Moreover, the initial states of the dynamical simulations are stationary solutions, provided by the selfconsistent procedure previously outlined. Figure 1a represents diagrammatically the π-electron force field, d n , in the absence of any disturbance by colored arrows. One may see that sites with the same parity possess arrows of equal size, meaning they experience the same potential. Fig. 1(b) shows this field for a polaron. As can be seen, sites far away from the quasiparticle continue to feel a uniform parity-based alternating force field. However, those near it have a local distortion on the force field. Figure 1c shows actual calculations of n of a system bearing a polaron. The conventional approach yields the red line. The effective method gives the blue one. The quasiparticle's presence breaks the uniformity of the potential leading to a local distortion on the force field. The inset presents n at several time snapshots. The undulation on the left of the red line is due to a breather. It is dependent on how the polaron is accelerated. The breather mode is discernible around the initial position of the polaron.

Results
The effective field is obtained by fitting n = d n + p n to a numerically calculated system with a polaron. We found that the polaron force field distortion, p n , is given by, in which, a 1 = −0.52 and a 2 = 0. 27 . Here we set v = a 3 A(t) with a 3 = 2.7 * 10 −3 . These parameters are robust and characterize polarons in different configurations. In practice, our methodology to estimate n begins by explicitly simulating polaron dynamics via the usual approach. Then, with this reference result, we use Eq. 10 to find B n , which translates into n through Eq. 11. Having this reference, the π force field is determined through the fitting of the parameters in Eq. 12. The fitted expression represents the general electronic contribution to the lattice dynamics of a conjugated polymer hosting a drifting polaron driven by an external field.
We present in Fig. 2 the results for the system's lattice and electronic pictures for comparison. Fig. 2a and c depict the lattice order parameter and charge density obtained through the conventional formalism. On the other hand, Fig. 2b and d refer to the same variables calculated using the effective potential.
The quasiparticle's trajectory consists of two phases: an initial acceleration and a regime of constant velocity. The first time interval ranges from 0 to 50 fs. The profile of this period is a direct outcome of the active presence of the electric field. The polaron begins to emit phonons in both directions. After 50 fs, the electric field goes off. Right after the phonon scattering event, the quasiparticles drift with a constant velocity until about 280 fs. At this time, the distortions reach the polaron initial position due to the periodic boundary condition. As soon as the quasiparticles return to this region, the polarons of both methodologies suffer from phonon collisions since they encounter breathers excited earlier. As the two structures collide, the signature distortion of the quasiparticle is momentarily muddled 7,17 .
(10) B n ≡ e −iγ A(t) k,s ′ ψ * k,s (n + 1, t)ψ k,s (n, t),  www.nature.com/scientificreports/ The trajectories delivered by both methodologies represent physically equivalent phenomena. In both cases, the polaron drifts by two distinct regimes, covering similar distances at the same time range. Additionally, the quasiparticles react identically when the acceleration ceases its influence. Figure 2c and d show the charge density time evolution. As expected, the motion of the charged region is the same as the local deformation that appears in the lattice order parameter. As in Fig. 2a-c, at the same time, the scattered phonons collide with the polaron. These events cause a disturbance in the localized lattice distortion. Because the quasiparticle is a coupled structure of lattice distortion and charge accumulation, a similar effect occurs in (c) and (d). However, since the scattered phonons are less energetic in (d), the disturbance in the charge density profile is smaller.
There is significant accordance regarding the energy dynamics. At the start of the simulations, until about 50 fs, the electronic energy rapidly increases with the applied electric field. Next, the curve ceases to grow and begins to follow an oscillatory trajectory until the end of the simulation. In addition, one may note that polaron moves slightly faster in the traditional simulation. This occurs because, in comparison with the usual approach, the π-field formalism generates a quasiparticle with more definite lattice deformation, increasing the carrier's transport inertia. For this reason, under same electric field conditions, the quasiparticle simulated via the usual approach will cover a slightly greater distance.
The equivalence between the two approaches also becomes apparent when the electronic band structure is analyzed. Figure 3 displays the time dependence of the system's electronic band around the gap with a polaron state. Here, we present the energy states from the two independent methods. The black-colored lines refer to the energy levels obtained with the standard simulation. Alternatively, the red-colored states result from the effective potential dynamics. In both cases, the lines represent electronic states. As can be seen, the methodologies deliver equivalent electronic descriptions of the system.
The overall behavior of the bands can be summarized in two distinct regimes. The first one starts at the beginning of the simulation and goes to about 50 fs. During this period, the energy levels suffer a short but continuous shrinkage towards the center of the gap. That is a direct consequence of the Peierls instability, which is triggered by the external field. The other regime initiates right after that and goes until the end of the simulation. The energy states display small fluctuations around an average value. This profile is a result of lattice disturbances caused by the polaron's continuous drift and the motion of scattered phonons 22 .
Besides the similarities, the band from the effective potential approach has more intense fluctuations. It is a consequence of the successive collisions with the normal oscillation mode visible in Fig. 2b. The energy bandgap of the conventional approach is found to be 1.81 ∓ 0.06 eV, while the effective field method gives 1.83 ∓ 0.1 eV.
The computation cost of the proposed model is a central point. As discussed in the methods, the approach eliminates the need to diagonalize the electronic Hamiltonian at each time step. Clearly, such optimization reduces the computational cost. All simulations were performed using the same machine. The results present a great drop in the execution time, averaging two orders of magnitude. Although such an estimate depends on the machine used, we expect no changes in the conclusion: the effective potential is two orders of magnitude faster than conventional simulations.

Conclusion
In conclusion, we introduced an alternative model to simulate π-conjugated systems. The approach describes the electronic contribution to lattice dynamics in terms of effective π-fields. We showed that in the absence of lattice disturbances the effective force field is uniform, behaving similarly to an alternating constant gravitational field. By extending the study to excited states, we showed that the presence of quasiparticles distorts the field locally. We ran extensive tests comparing it with the usual formalism to estimate the model's reliability. The comparison between the modelings took place with dynamical simulations of the drifting of polaron under the influence of external electric fields. The two methodologies agree on the polaron's trajectory, the electronic www.nature.com/scientificreports/ band structure, and the phonon scattering. In that sense, the adoption of this new strategy causes no loss in the description's physical accuracy. In addition, we also tested the model's limits by enforcing different electric field regimes. Even under such conditions, the model coincides tightly with reference simulations. We also report a sharp time drop of two orders of magnitude compared with the standard approach. That is due to freeing the calculations from the constraint of diagonalization at each time in the evolution process. Our approach is an alternative to simulate much larger π-systems, where π-electron interactions are considered, and dynamics of multiple polarons become possible.

Data availability
The datasets used and/or analysed during the current study available from the corresponding author on reasonable request.